main_code.utf8.md

if(!require(ggplot2))install.packages("ggplot2"); library(ggplot2)

## Loading required package: ggplot2

if(!require(ggmap)) install.packages("ggmap"); library(ggmap)

## Loading required package: ggmap

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

# ggmap(get_map(location='south korea', zoom=7))
#error 발생 : google maps API terms of Service가 변경됨

if(!requireNamespace("devtools")) install.packages("devtools")

## Loading required namespace: devtools

devtools::install_github("dkahle/ggmap", ref = "tidyup")

## Skipping install of 'ggmap' from a github remote, the SHA1 (2d756e5e) has not changed since last install.
##   Use `force = TRUE` to force installation

# 좌, 우, 위, 아래의 값을 지정해 주어야 지도 그림이 그려짐
kor <- c(left = 124, bottom = 33, right = 132, top = 39)
map <- get_stamenmap(kor, zoom = 6, maptype = "toner-lite")

## Source : http://tile.stamen.com/toner-lite/6/54/24.png

## Source : http://tile.stamen.com/toner-lite/6/55/24.png

## Source : http://tile.stamen.com/toner-lite/6/54/25.png

## Source : http://tile.stamen.com/toner-lite/6/55/25.png

ggmap(map)

map <- get_stamenmap(kor, 
                     zoom = 6, maptype = "watercolor")

## Source : http://tile.stamen.com/watercolor/6/54/24.jpg

## Source : http://tile.stamen.com/watercolor/6/55/24.jpg

## Source : http://tile.stamen.com/watercolor/6/54/25.jpg

## Source : http://tile.stamen.com/watercolor/6/55/25.jpg

ggmap(map)

if(!require(dplyr)) install.packages("dplyr");library(dplyr)

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

if(!require(forcats)) install.packages("forcats");library(forcats)

## Loading required package: forcats

setwd("C:/Users/USER/Dropbox/R-work/vis/day2/extra")

load(file='crime.rda')
head(crime)

##                      time     date hour premise            offense  beat
## 82729 2010-01-01 15:00:00 1/1/2010    0     18A             murder 15E30
## 82730 2010-01-01 15:00:00 1/1/2010    0     13R            robbery 13D10
## 82731 2010-01-01 15:00:00 1/1/2010    0     20R aggravated assault 16E20
## 82732 2010-01-01 15:00:00 1/1/2010    0     20R aggravated assault  2A30
## 82733 2010-01-01 15:00:00 1/1/2010    0     20A aggravated assault 14D20
## 82734 2010-01-01 15:00:00 1/1/2010    0     20R           burglary 18F60
##           block    street type suffix number   month    day
## 82729 9600-9699   marlive   ln      -      1 january friday
## 82730 4700-4799 telephone   rd      -      1 january friday
## 82731 5000-5099  wickview   ln      -      1 january friday
## 82732 1000-1099   ashland   st      -      1 january friday
## 82733 8300-8399    canyon           -      1 january friday
## 82734 9300-9399     rowan   ln      -      1 january friday
##                       location           address       lon      lat
## 82729    apartment parking lot   9650 marlive ln -95.43739 29.67790
## 82730 road / street / sidewalk 4750 telephone rd -95.29888 29.69171
## 82731        residence / house  5050 wickview ln -95.45586 29.59922
## 82732        residence / house   1050 ashland st -95.40334 29.79024
## 82733                apartment       8350 canyon -95.37791 29.67063
## 82734        residence / house     9350 rowan ln -95.54830 29.70223

table(crime$offense)

## 
## aggravated assault         auto theft           burglary 
##               7177               7946              17802 
##             murder               rape            robbery 
##                157                378               6298 
##              theft 
##              46556

# define helper
`%notin%` <- function(lhs, rhs) !(lhs %in% rhs)

# reduce crime to violent crimes in downtown houston
violent_crimes <- crime %>% 
  filter(
    offense %notin% c("auto theft", "theft", "burglary"),
    -95.39681 <= lon & lon <= -95.34188,
    29.73631 <= lat & lat <=  29.78400
  ) %>% 
  mutate(
    offense = fct_drop(offense), # refactor 
    offense = fct_relevel(offense, 
                          c("robbery", "aggravated assault", "rape", "murder") #factor order
    )
  )

head(violent_crimes)

##                  time     date hour premise            offense  beat
## 1 2010-01-01 20:00:00 1/1/2010    5     20R aggravated assault  2A10
## 2 2010-01-01 21:00:00 1/1/2010    6     20R               rape 10H10
## 3 2010-01-03 09:00:00 1/2/2010   18     18T            robbery  1A20
## 4 2010-01-04 05:00:00 1/3/2010   14     210            robbery  6B10
## 5 2010-01-05 19:00:00 1/5/2010    4     20A aggravated assault  2A30
## 6 2010-01-06 05:00:00 1/5/2010   14     070 aggravated assault 10H40
##       block    street    type suffix number   month      day
## 1 2000-2099     keene              -      1 january   friday
## 2 2300-2399   sperber      ln      -      1 january   friday
## 3   400-499      gray              W      1 january saturday
## 4 5500-5599     north fwy ser      E      1 january   sunday
## 5 2200-2299 white oak      dr      -      1 january  tuesday
## 6 1000-1099   alabama      st      -      1 january  tuesday
##                            location            address       lon      lat
## 1                 residence / house         2050 keene -95.36348 29.77802
## 2                 residence / house    2350 sperber ln -95.34817 29.75505
## 3 strip business center parking lot           450 gray -95.37596 29.75143
## 4            restaurant / cafeteria 5550 north fwy ser -95.36327 29.76328
## 5                         apartment  2250 white oak dr -95.38217 29.78002
## 6                 convenience store    1050 alabama st -95.37930 29.73745

table(violent_crimes$offense)

## 
##            robbery aggravated assault               rape 
##                312                362                 22 
##             murder 
##                  5

# use qmplot to make a scatterplot on a map
qmplot(lon, lat, data = violent_crimes, maptype = "toner-lite", color = I("red"))

## Using zoom = 14...

## Source : http://tile.stamen.com/terrain/14/3850/6770.png

## Source : http://tile.stamen.com/terrain/14/3851/6770.png

## Source : http://tile.stamen.com/terrain/14/3852/6770.png

## Source : http://tile.stamen.com/terrain/14/3853/6770.png

## Source : http://tile.stamen.com/terrain/14/3850/6771.png

## Source : http://tile.stamen.com/terrain/14/3851/6771.png

## Source : http://tile.stamen.com/terrain/14/3852/6771.png

## Source : http://tile.stamen.com/terrain/14/3853/6771.png

## Source : http://tile.stamen.com/terrain/14/3850/6772.png

## Source : http://tile.stamen.com/terrain/14/3851/6772.png

## Source : http://tile.stamen.com/terrain/14/3852/6772.png

## Source : http://tile.stamen.com/terrain/14/3853/6772.png

## Source : http://tile.stamen.com/terrain/14/3850/6773.png

## Source : http://tile.stamen.com/terrain/14/3851/6773.png

## Source : http://tile.stamen.com/terrain/14/3852/6773.png

## Source : http://tile.stamen.com/terrain/14/3853/6773.png

qmplot(lon, lat, data = violent_crimes, maptype = "toner-lite", geom = "density2d", color = I("red"))

## Using zoom = 14...

robberies <- violent_crimes %>% filter(offense == "robbery")

qmplot(lon, lat, data = violent_crimes, geom = "blank", 
       zoom = 15, maptype = "toner-background", darken = .7, legend = "topleft") +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = .3, color = NA)+
  scale_fill_gradient2("Robbery\nPropensity", low = "white", mid = "yellow", high = "red", midpoint = 650)

## 49 tiles needed, this may take a while (try a smaller zoom).

## Source : http://tile.stamen.com/terrain/15/7700/13541.png

## Source : http://tile.stamen.com/terrain/15/7701/13541.png

## Source : http://tile.stamen.com/terrain/15/7702/13541.png

## Source : http://tile.stamen.com/terrain/15/7703/13541.png

## Source : http://tile.stamen.com/terrain/15/7704/13541.png

## Source : http://tile.stamen.com/terrain/15/7705/13541.png

## Source : http://tile.stamen.com/terrain/15/7706/13541.png

## Source : http://tile.stamen.com/terrain/15/7700/13542.png

## Source : http://tile.stamen.com/terrain/15/7701/13542.png

## Source : http://tile.stamen.com/terrain/15/7702/13542.png

## Source : http://tile.stamen.com/terrain/15/7703/13542.png

## Source : http://tile.stamen.com/terrain/15/7704/13542.png

## Source : http://tile.stamen.com/terrain/15/7705/13542.png

## Source : http://tile.stamen.com/terrain/15/7706/13542.png

## Source : http://tile.stamen.com/terrain/15/7700/13543.png

## Source : http://tile.stamen.com/terrain/15/7701/13543.png

## Source : http://tile.stamen.com/terrain/15/7702/13543.png

## Source : http://tile.stamen.com/terrain/15/7703/13543.png

## Source : http://tile.stamen.com/terrain/15/7704/13543.png

## Source : http://tile.stamen.com/terrain/15/7705/13543.png

## Source : http://tile.stamen.com/terrain/15/7706/13543.png

## Source : http://tile.stamen.com/terrain/15/7700/13544.png

## Source : http://tile.stamen.com/terrain/15/7701/13544.png

## Source : http://tile.stamen.com/terrain/15/7702/13544.png

## Source : http://tile.stamen.com/terrain/15/7703/13544.png

## Source : http://tile.stamen.com/terrain/15/7704/13544.png

## Source : http://tile.stamen.com/terrain/15/7705/13544.png

## Source : http://tile.stamen.com/terrain/15/7706/13544.png

## Source : http://tile.stamen.com/terrain/15/7700/13545.png

## Source : http://tile.stamen.com/terrain/15/7701/13545.png

## Source : http://tile.stamen.com/terrain/15/7702/13545.png

## Source : http://tile.stamen.com/terrain/15/7703/13545.png

## Source : http://tile.stamen.com/terrain/15/7704/13545.png

## Source : http://tile.stamen.com/terrain/15/7705/13545.png

## Source : http://tile.stamen.com/terrain/15/7706/13545.png

## Source : http://tile.stamen.com/terrain/15/7700/13546.png

## Source : http://tile.stamen.com/terrain/15/7701/13546.png

## Source : http://tile.stamen.com/terrain/15/7702/13546.png

## Source : http://tile.stamen.com/terrain/15/7703/13546.png

## Source : http://tile.stamen.com/terrain/15/7704/13546.png

## Source : http://tile.stamen.com/terrain/15/7705/13546.png

## Source : http://tile.stamen.com/terrain/15/7706/13546.png

## Source : http://tile.stamen.com/terrain/15/7700/13547.png

## Source : http://tile.stamen.com/terrain/15/7701/13547.png

## Source : http://tile.stamen.com/terrain/15/7702/13547.png

## Source : http://tile.stamen.com/terrain/15/7703/13547.png

## Source : http://tile.stamen.com/terrain/15/7704/13547.png

## Source : http://tile.stamen.com/terrain/15/7705/13547.png

## Source : http://tile.stamen.com/terrain/15/7706/13547.png

qmplot(lon, lat, data = violent_crimes, maptype = "toner-background", color = offense) + 
  facet_wrap(~ offense)

## Using zoom = 14...

#유럽지도 가져오기
europe <- c(left = -12, bottom = 35, right = 30, top = 63)
get_stamenmap(europe, zoom = 5) %>% ggmap()

## Source : http://tile.stamen.com/terrain/5/14/8.png

## Source : http://tile.stamen.com/terrain/5/15/8.png

## Source : http://tile.stamen.com/terrain/5/16/8.png

## Source : http://tile.stamen.com/terrain/5/17/8.png

## Source : http://tile.stamen.com/terrain/5/18/8.png

## Source : http://tile.stamen.com/terrain/5/14/9.png

## Source : http://tile.stamen.com/terrain/5/15/9.png

## Source : http://tile.stamen.com/terrain/5/16/9.png

## Source : http://tile.stamen.com/terrain/5/17/9.png

## Source : http://tile.stamen.com/terrain/5/18/9.png

## Source : http://tile.stamen.com/terrain/5/14/10.png

## Source : http://tile.stamen.com/terrain/5/15/10.png

## Source : http://tile.stamen.com/terrain/5/16/10.png

## Source : http://tile.stamen.com/terrain/5/17/10.png

## Source : http://tile.stamen.com/terrain/5/18/10.png

## Source : http://tile.stamen.com/terrain/5/14/11.png

## Source : http://tile.stamen.com/terrain/5/15/11.png

## Source : http://tile.stamen.com/terrain/5/16/11.png

## Source : http://tile.stamen.com/terrain/5/17/11.png

## Source : http://tile.stamen.com/terrain/5/18/11.png

## Source : http://tile.stamen.com/terrain/5/14/12.png

## Source : http://tile.stamen.com/terrain/5/15/12.png

## Source : http://tile.stamen.com/terrain/5/16/12.png

## Source : http://tile.stamen.com/terrain/5/17/12.png

## Source : http://tile.stamen.com/terrain/5/18/12.png

kor <- c(left = 124, bottom = 33, right = 132, top = 39)
map <- get_stamenmap(kor, zoom = 7, maptype = "toner-lite")

## Source : http://tile.stamen.com/toner-lite/7/108/48.png

## Source : http://tile.stamen.com/toner-lite/7/109/48.png

## Source : http://tile.stamen.com/toner-lite/7/110/48.png

## Source : http://tile.stamen.com/toner-lite/7/108/49.png

## Source : http://tile.stamen.com/toner-lite/7/109/49.png

## Source : http://tile.stamen.com/toner-lite/7/110/49.png

## Source : http://tile.stamen.com/toner-lite/7/108/50.png

## Source : http://tile.stamen.com/toner-lite/7/109/50.png

## Source : http://tile.stamen.com/toner-lite/7/110/50.png

## Source : http://tile.stamen.com/toner-lite/7/108/51.png

## Source : http://tile.stamen.com/toner-lite/7/109/51.png

## Source : http://tile.stamen.com/toner-lite/7/110/51.png

plot(map)

wifi <- read.csv('wifi.csv', header=T) #wifi.csv
head(wifi)

##   company      lat      lon
## 1      KT 37.74417 128.9056
## 2      KT 37.72806 128.9543
## 3      KT 37.75710 128.8900
## 4      KT 37.74769 128.8840
## 5      KT 37.74866 128.9073
## 6      KT 37.74281 128.8827

ggmap(map) + geom_point(data=wifi, aes(x=lon, y=lat, color=company))

ggmap(map) + stat_density_2d(data=wifi, aes(x=lon, y=lat))

ggmap(map) + 
  stat_density_2d(data=wifi, aes(x=lon, y=lat, fill=..level.., alpha=..level..), geom='polygon', size=2, bins=30)

#일단 aes 안에 들어 있는 fill은 문자 그대로 색깔로 채우라는 뜻입니다.
# ..level..은 레벨(level)이 높을수록,  예를 들어 기압이 높거나 고도가 높을수록
# 더 진한 색깔을 칠하라는 뜻입니다.
# alpha는 투명도를 나타냅니다.
# 역시 같은 원리로 레벨이 높으면 불투명하게(색이 더 잘 드러나게)
# 칠하고 낮을 때는 투명하게(희미하게) 칠하라는 의미입니다

# geom='polygon'에서 polygon은 다각형이라는 뜻입니다.
# 기본은 위에서 본 것처럼 선(線)으로 돼 있는데 그것 말고 도형으로 그리라고 명령을 준 겁니다.
# 선을 기준으로 size는 선 굵기, bins는 선 간격이라는 뜻입니다.
# 이름이 이렇게 붙은 건 원래 점을 이은 선을 기준으로 삼고 있는 까닭입니다.

p <- ggmap(map) + stat_density_2d(data=wifi, aes(x=lon, y=lat, fill=..level.., alpha=..level..), geom='polygon', size=7, bins=28)
p

p<-p + scale_fill_gradient(low='yellow', high='red', guide=F)
p

p<-p + scale_alpha(range=c(0.02, 0.8), guide=F)
p

airport <- read.csv('airport.csv', header=T) #airport.csv

route <- read.csv('route.csv', header=T) #route.csv

head(airport)

##   airport iata     lon     lat
## 1    강릉  KAG 128.944 37.7536
## 2    광주  KWJ 126.809 35.1264
## 3    군산  KUV 126.616 35.9038
## 4    김포  GMP 126.791 37.5583
## 5    대구  TAE 128.659 35.8941
## 6    목포  MPK 126.380 34.7589

head(route)

##   id airport     lon     lat
## 1  1     CJJ 127.499 36.7166
## 2  7     CJJ 127.499 36.7166
## 3 45     CJJ 127.499 36.7166
## 4 77     CJJ 127.499 36.7166
## 5  2     CJJ 127.499 36.7166
## 6  8     CJJ 127.499 36.7166

ggmap(map) + geom_point(data=airport, aes(x=lon, y=lat, size=3))

p <- ggmap(map) + geom_point(data=airport, aes(x=lon, y=lat, size=3))

p<-p + geom_line(data=route, aes(x=lon, y=lat, group=id))
p

CJU_route<-route[route$airport=="CJU","id"]
route1 <- route[route$id %in% CJU_route,]
p + geom_line(data=route1, aes(x=lon, y=lat, group=id, size=2), col='gold')

# Q. 기상관측소 데이터를 찾아서 지도에 그려보세요.
# stn<- read.csv(choose.files())

if(!require(raster))install.packages("raster"); library(raster)

## Loading required package: raster

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

korea <- getData('GADM', country='kor', level=2)
plot(korea)

ggplot(data=korea, aes(x=long, y=lat, group=group)) +
  geom_polygon(fill='white', color='black')

## Regions defined for each Polygons

# 잘 보시면 예전에 못 보던 게 하나 보입니다.
# group : 선을 모아 면을 만들 때도 짝을 맞춰주는 속성
# 이 지도는 옛날 행정구역임

if(!require(rgdal))install.packages("rgdal"); library(rgdal)

## Loading required package: rgdal

## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/USER/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/USER/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1

if(!require(maptools))install.packages("maptools"); library(maptools)

## Loading required package: maptools

## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()

korea <- readOGR('SIG_201703/TL_SCCO_SIG.shp')#'TL_SCCO_SIG.shp'

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\USER\Dropbox\R-work\vis\day2\extra\SIG_201703\TL_SCCO_SIG.shp", layer: "TL_SCCO_SIG"
## with 250 features
## It has 3 fields

ggplot() + geom_polygon(data=korea, aes(x=long, y=lat, group=group), fill='white', color='black')

## Regions defined for each Polygons

#19대 대통령 선거 시군구별 득표율
result <- read.csv('result.csv', header=T) #'result.csv'
head(result)

##   Sido    Sigun    id  moon  hong   ahn   yoo  shim   etc  X1st X1st_C
## 1 서울   종로구 11110 0.416 0.218 0.218 0.073 0.070 0.005 0.416   Moon
## 2 서울     중구 11140 0.412 0.217 0.235 0.071 0.060 0.005 0.412   Moon
## 3 서울   용산구 11170 0.393 0.239 0.217 0.080 0.066 0.004 0.393   Moon
## 4 서울   성동구 11200 0.428 0.200 0.226 0.078 0.064 0.004 0.428   Moon
## 5 서울   광진구 11215 0.441 0.194 0.221 0.072 0.069 0.004 0.441   Moon
## 6 서울 동대문구 11230 0.421 0.219 0.227 0.064 0.064 0.005 0.421   Moon
##    X2nd X2nd_C   gap final
## 1 0.218   Hong 0.198 0.416
## 2 0.235    Ahn 0.178 0.412
## 3 0.239   Hong 0.155 0.393
## 4 0.226    Ahn 0.203 0.428
## 5 0.221    Ahn 0.220 0.441
## 6 0.227    Ahn 0.194 0.421

head(korea)

##   SIG_CD    SIG_ENG_NM SIG_KOR_NM
## 0  11110     Jongno-gu     종로구
## 1  11140       Jung-gu       중구
## 2  11170    Yongsan-gu     용산구
## 3  11200  Seongdong-gu     성동구
## 4  11215   Gwangjin-gu     광진구
## 5  11230 Dongdaemun-gu   동대문구

korea1<-merge(korea, result, by.x="SIG_CD", by.y="id")
head(korea1)

##   SIG_CD    SIG_ENG_NM SIG_KOR_NM Sido    Sigun  moon  hong   ahn   yoo
## 1  11110     Jongno-gu     종로구 서울   종로구 0.416 0.218 0.218 0.073
## 2  11140       Jung-gu       중구 서울     중구 0.412 0.217 0.235 0.071
## 3  11170    Yongsan-gu     용산구 서울   용산구 0.393 0.239 0.217 0.080
## 4  11200  Seongdong-gu     성동구 서울   성동구 0.428 0.200 0.226 0.078
## 5  11215   Gwangjin-gu     광진구 서울   광진구 0.441 0.194 0.221 0.072
## 6  11230 Dongdaemun-gu   동대문구 서울 동대문구 0.421 0.219 0.227 0.064
##    shim   etc  X1st X1st_C  X2nd X2nd_C   gap final
## 1 0.070 0.005 0.416   Moon 0.218   Hong 0.198 0.416
## 2 0.060 0.005 0.412   Moon 0.235    Ahn 0.178 0.412
## 3 0.066 0.004 0.393   Moon 0.239   Hong 0.155 0.393
## 4 0.064 0.004 0.428   Moon 0.226    Ahn 0.203 0.428
## 5 0.069 0.004 0.441   Moon 0.221    Ahn 0.220 0.441
## 6 0.064 0.005 0.421   Moon 0.227    Ahn 0.194 0.421

SIG_CD<-paste(unique(korea1$SIG_CD))

f.korea <- fortify(korea) #maptools package

## Regions defined for each Polygons

f.korea$SIG_CD <-NA

for (i in 1:250){
  f.korea[f.korea$id==(i-1),"SIG_CD"]=SIG_CD[i]  
}

korea2 <- merge(f.korea, result, by.x="SIG_CD", by.y="id")

p<-ggplot() + geom_polygon(data=korea2, aes(x=long, y=lat, group=group, fill=moon))
p + scale_fill_gradient(low='white', high='#004ea2')

if(!require(viridis))install.packages("viridis");library(viridis)

## Loading required package: viridis

## Loading required package: viridisLite

p<-ggplot() + geom_polygon(data=korea2, aes(x=long, y=lat, group=group, fill=moon))
p + scale_fill_viridis()

# 기존에 있는 파레트를 이용해서 그려보자.

p + scale_fill_viridis(direction=-1) + theme_void() + guides(fill=F)

# 색을 반대로, 회색을 지우고, 범례도 지워보자

# 다른 후보들의 그림을 그려보세요.
# 문재인 후보와 홍준표 후보의 차이를 시각적으로 살펴보세요.

main_code.R

USER

2019-06-24